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Abstract 

Diffusion Monte Carlo results for the phonon-roton excitation branch in bulk 
liquid ^He at zero temperature are presented. The sign problem associated 
to the excited wave function has been dealt both with the fixed-node ap- 
proximation and the released-node technique. The upper bounds provided 
by the fixed-node approximation are shown to become exact when using the 
released-node method. An excellent agreement with experimental data is 

achieved both at the equilibrium and near the freezing densities. 
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The physical nature of the excitations in superfluid ^He at low momenta is still nowadays 
not completely well understood [|I|, in contrast with the vast knowledge of its static prop- 
erties. In a pioneering work, Landau [0 proposed a model for the elementary excitations 
to explain the superfluidity of liquid ^He. From a different point of view, Bogoliubov 
calculated the excitation spectrum of a weakly interacting Bose gas where the condensate 
fraction, i.e., the fraction of particles in the zero momentum state, plays an explicit role. 
Both theories predict a continuous dispersion curve which starts with a phonon excitation, 
reaches a first maximum (maxon), lowers to a local minimum (roton), and then grows up as 
the energy of a free particle. From general ideas on the nature of the excitations in an inter- 
acting Bose fiuid, Feynman proposed the first microscopic approach to the problem. The 
Feynman trial wave function provides a qualitative description of the excitation spectrum 
but fails in reproducing the roton energy by a factor two. Later on, Feynman and Cohen 
included backfiow correlations in the trial wave function, reducing in one half the differences 
between the experimental data and the original Feynman's prediction. Following Feynman's 
language, the phonon-roton branch corresponds to collective density-like excitations where 
the condensate fraction does not enter in an explicit way. Recently, it has been argued by 
Clyde and Criffin ^ that the continuous spectrum results from the superposition of den- 
sity excitations, dominating in the phonon region, and single-particle excitations, important 
in the roton minimum. This theory has emerged after the experimental determination of 
the temperature dependence of the excitation spectrum which shows the phonon peak in 
the dynamic structure function S{q, u) at both sides of the A-transition whereas the roton 
practically disappears in the normal phase. 

In the last years a considerable effort has been made to improve quantitatively the 
microscopic predictions for the phonon-roton excitation spectrum e{q). Manousakis and 
Pandharipande ^ calculated e{q) by means of the correlated basis function (CBF) method 
using a basis of Feynman-Cohen states. The variational Monte Carlo (VMC) method using 
shadow wave functions has also proved to be quantitatively quite efficient in the calculation 
of e{q) l^l in spite of its approximate description of the ground state. The application of ab 
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initio Monte Carlo methods to this problem has been, however, severely hindered by the sign 
problem associated to the excited wave function. Only recently, Boninsegni and Ceperley 
have calculated 6{q) by means of a path integral Monte Carlo (PIMC) calculation of S{q,uj) 
from a Laplace inversion of the imaginary-time correlation factor S{q, t). However, for noisy 
data this inversion is an ill-posed problem that prevents a model-independent determination. 

In the present work, we present a zero temperature calculation of the phonon-roton 
spectrum using the diffusion Monte Carlo (DMC) method. In this method, the imaginary- 
time Schrodinger equation for the function /(R, t) = '?/;t(R)$(R, t), 

- = -DVifiR, t) + DVr(F(R)/(R, t)) + (E4R) - E)f{R, t) , (1) 

is solved stochastically, $(R, t) and iPtCR) being the wave function of the system and a trial 
wave function used for importance sampling, respectively. In the above equation, -E'l(R) = 
ipT'^i'R)H'ipT(R) is the local energy and F(R) = 2'?/'^"'^(R)Vr'?/'t(R) acts as a quantum drift 
force; D = fi^ /2m, with m the mass of the particles, and R stands for the 3A^-coordinate 
vector of the particles. In the asymptotic regime, /(R, t oo) — > ipxCR-) ^(R) where 
$(R) corresponds to the lowest energy eigenstate of the system not orthogonal to the trial 
wave function i/jtCR). For a thorough description of the DMC method see, for instance. 
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The wave function ipCR.) corresponding to a phonon-roton excitation is an eigenstate of 
the momentum operator. As pointed out by Feynman [Q, this requirement is achieved with 
the simple model wave function 

N 

V^^(R) = ^e^''->^°(R) , (2) 

i=l 

tp^(R) being the ground-state wave function. The first correction to iIj^(R), originally 
proposed by Feynman and Cohen ||^, includes backfiow correlations. In this case, the excited 
wave function is given by 

TV 

r^{R) = J2e'''-'^ rm , (3) 

i=l 



where 



= ri + ^r]{rij) . (4) 

The inclusion of backflow correlations improves appreciably the variational results of e{q) 
with respect to the Feynman's choice but the quantitative agreement with experiment is 
still poor, specially near the roton minimum. 

In the DMC implementation a real probability distribution function /(R, t) is suitable. 
Therefore, we choose as importance sampling wave function ipxi^) the superposition of two 
excitations of momenta q and — q which are degenerate in energy: 

N 

^^(R) = ^cos(q-r,)^°(R) (5) 



at Feynman's level, and 



1=1 



N 



^|^(R) = ^cos(q-i-,)^^(R) (6) 



i=l 



when backflow correlations are included. 

In a first step, the sign problem associated to the excited wave function has been dealt 
within the framework of the fixed-node (FN) approximation [^]. The FN approximation 
provides an upper bound to the exact value and has been extensively used in the study 
of fermionic systems. In what concerns the calculation of excited states, we have recently 
used the FN-DMC method to study a vortex excitation in two-dimensional superfluid ^He 
1^ . Within FN, the sign problem is avoided by imposing the nodal surface of the trial 



wave function to the excited Monte Carlo wave function. The problem is hence mapped 
onto a bosonic calculation with a well defined density probability function. Beyond this 
approximation, one can remove the nodal constraint of the FN approach using a released- 
node technique (RN) |T3[. This method makes use of an auxiliary guiding function ipgili), 
positively defined everywhere. Random walkers are then allowed to cross the nodal surface 
(becoming negative) and to survive for a finite lifetime tr- The RN approach provides 
the exact eigenvalue in the limit of large surviving times, but presents the disadvantage 



of becoming numerically unstable in the limit tr ^ oo where the number of positive and 
negative walkers become similar. 

The success of the method, i.e., the achievement of an asymptotic regime before the 
growth of the statistical errors, is closely related to the quality of iPt(R) and ipg{H). The 
guiding function ipgill) has to approach |'^t(R)| away from the nodal surface and must be 
non-zero in the nodes to make possible the flux of walkers through it. We have taken the 
simple model 

V',(R) = (^T(R)' + a')'^' , (7) 

which satisfies both requirements for a proper choice of the value of the parameter a. In 
fact, the value given to the parameter a, which has to be of the same order of magnitude 
than the mean value of |'?/'t(R')|, governs the flux of walkers through the nodes. Therefore, 
the relaxation time of the release process is a function of a: when the value of a is increased 
the relaxation time is reduced and vice versa. 

The released-node energy is obtained projecting out the excited state modeled by ipri^)- 
This projection is carried out assigning to each walker a weight iy(R) given by 

(j(R) being +1 (—1) for an even (odd) number of crossings. The released-node energy is 
thus determined through 

where the sums are extended to all the surviving walkers with a lifetime less than tr, and 

In our computation we have considered a simulation box with 108 ^He atoms interacting 
through the IIFD-B(IIE) Aziz potential [0], and at two densities, po = 0.365 cr~^ and 
Pp = 0.438 o"^^ (cr = 2.556 A). The density po corresponds to the equilibrium density, and 
Pp is close to the freezing density. The ground-state correlations have been modeled by a 
two-body wave function originally proposed by Reatto [|T^ , 



(10) 



which we have previously employed in ground-state calculations |Tl]|. The function ri{r), 



entering in the backflow wave function (Hf), has been chosen to be a gaussian 



r]{r) = A exp 



^ 9 

r — Th ^ 



(11) 
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as in the variational calculation of Ref. 0. The values of the parameters in Eqs. (|lOip!lD 



are b = 1.20 a, L = 0.2, A = 0.6 (x, A = 2.0 a, A = 0.15, n = 0.8 a, and Ub = 0.44 a, which 
are roughly the optimal ones at the equilibrium density po- The parameter a appearing in 
the guiding wave function ipgC^) (0) been taken as a = 3.0 for all the q values, and the 
largest lifetime tr used corresponds to 190 DMC sweeps. The same set of parameters have 
been used at the highest density pp. 

As mentioned before, the released-node energy estimation would be exact in the limit of 
large lifetimes. However, the computational effort to simultaneously enlarge tr and main- 
tain the statistical fluctuations into an acceptable level grows with tr- Having estimated a 
reasonable upper limit of t^, given the present computational resources, we can study the 
influence of the excited trial wave function in the RN energy. As a general trend, if the RN 
energy does not reach a constant regime within tr, an improved model for the excited trial 
wave function should be used. The empirical way in which we have studied the asymptotic 
regime of the RN energy is by fitting the function 

E{tr) = E^ + C e"*-/" (12) 

to the largest tr values. In the interpretation of the MC results, we have followed the guide- 
line of accepting only the RN values that do not present discrepancies between the largest 
tr data and the asymptotic limit Eao- The fit (|T^) has been used to decide whether to trust 
or not the MC values but not to provide the asymptotic limit. 

We have verified that at po, and for values q < 2.5 A^^, the RN energies using the 
Feynman wave function (|]) do reach the expected constant regime, the difference between 



the largest tr calculated energy and the value of £"00 predicted by the x^-fit (Q) being less 
than the statistical error. This is not the case for pp. At this high density, that agreement 
only subsists for the lowest q value and for the value of q nearest to the roton. For the other 
values of q we have had to include backfiow correlations (^) to reach the asymptotic limit. 
This fact is illustrated in Fig. 1, where the excitation energy 

^ iMqmm) _ {rT\H\^') 

per particle is plotted as a function of tr for q = 1.11 A^^ and q = 1.84 A~^. Near the roton, 
q = 1.84 A~^, both the Feynman and backfiow results show a coincident asymptotic value 
without a significant slope. At g = 1.11 A^^, near the maxon energy, the situation is clearly 
different: the backfiow results have reached a constant behaviour whereas the Feynman ones 
show a slow approach to the correct value, which has not been achieved yet for the maximum 
value of tr- The latter behaviour is also observed for the highest value of g (g ^ 2.6 A~^) 
both at po and pp. In this case, the inclusion of backfiow correlations in the wave function 
is not enough to eliminate the bias and a significant difference exists between the largest tr 
energy and the asymptotic value predicted by the numerical fit (|12|). 

The released-node mechanism suppresses the fixed-node constraints and drives the cal- 
culation to the exact excitation energy, as shown in Table I. In the table, fixed-node values 
using ip^ and ipp^ , and the released-node estimation are compared with experimental data 
||T6| at the equilibrium density po- The FN results with backfiow correlations improve the 
Feynman ones for the three values of g in a magnitude which depends on g. Thus, the inclu- 
sion of backfiow correlations seems slightly more relevant in the roton than in the maxon. 
On the other hand, the RN excitation energies agree with the experimental data for the 
three values of g within the statistical errors. 



In Fig. 2 the RN excitation energies are compared with the experimental spectrum ||T6 
at Po. The RN results correspond, for each g, to the last point in the release process, the 
error bars being only the statistical errors. As commented before, the systematic errors 
are less than the statistical ones except for the highest g result (g = 2.58 A~^). For this 



latter value of q we also report an estimation coming from the extrapolation supplied by 
the fit (|12|). Apart from this point, where the RN method shows the shortcomings of the 
backflow wave function at a so high value of g, the agreement between the RN results and 
the experiment is excellent. As a matter of comparison, the FN results using ip^^ are also 
plotted. It is worth noticing the difference between the FN energies in the maxon and in the 
roton regimes; the roton is reproduced quite accurately whereas in the maxon the backflow 
correlations overestimate appreciably the excitation energies. At the highest g, where the 
spectrum bends down, the FN energy is quite far from the experimental data. 

As is well known from neutron scattering data, the location and depth of the roton 
minimum depends on the density. Thus, when the density increases the roton appears shifted 
to higher momenta and its energy decreases. The energies in the maxon region increase with 
the density but in an amount not so well experimentally known as in the roton. In Fig. 3, we 



report the RN excitation energies at pp in comparison with experimental data IT^. There is, 
again, a good agreement between theory and experiment within the statistical errors except 
at the highest q evaluated {q = 2.74 A~^) where the RN energy has not reached a constant 
value inside the release interval. 

In conclusion, we have shown that the diffusion Monte Carlo method in conjunction 
with the fixed-node technique, and more specially, with the released-node method provides 
a very useful tool to study excitations in correlated quantum many body systems like liquid 
'^He. The results for e{q) are in an excellent quantitative agreement with experimental data, 
both at the equilibrium and near the freezing densities, improving previous variational and 
CBF results and providing an exact description of one of the oldest and hardest problems 
in the study of quantum fluids from a microscopical viewpoint. Possible applications of 
the RN-DMC method would be the ripplon excitations in a free liquid ^He surface or the 
determination of the excitation energy of a single impurity in bulk liquid ^He. On the other 
hand, the interpretation of the roton excitation as a single-particle mode deserves further 
theoretical work from a microscopical viewpoint |TB . 
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FIGURES 

FIG. 1. Excitation energies per particle as a function of the lifetime tr at pp. The full circles 
are obtained using ^/^^^ and the diamonds using ip^. 

FIG. 2. Phonon-roton spectrum at the equilibrium density pQ. The full circles are the RN 
results and the diamonds correspond to a FN calculation with . The open square, which has 
been slightly shifted to the right for clarity, is the result of the extrapolation with the fit (0). The 
solid line is the experimental data from Ref. [J16U. 

FIG. 3. Phonon-roton spectrum at the density pp. Same notation as in Fig. 2. The experi- 
mental data is from Ref. [^]. 
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TABLES 

TABLE L Excitation energies at po in comparison with experimental data. The FN-ip^ and 
FN-^y^ columns are the fixed- node energies using ipj^ and ^7?^, respectively. The RN column 



corresponds to the released- node estimation. Experimental data is taken from Ref. |16| 



FN-^|: (K) 



FN-^I?^ (K) 



RN (K) 



Expt. (K) 



0.369 
L106 
1.844 



7.56 ± 0.49 
18.47± 0.49 
13.82± 0.54 



7.24 ± 0.38 
16.52 ± 0.43 
10.37 ± 0.59 



7.02 ± 0.49 
13.82 ± 0.43 
9.18 ± 0.59 



7.0 
13.8 
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